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INTRODUCTION  TO  VOLUME  II 


This  volume  consists  of  a  series  of  appendices  whose  contents  could 
not  be  included  in  the  body  of  the  report  without  breaking  the  continuity. 
The  appendices  describe  topics  that  were  studied  during  the  course  of  the 
current  contract  but  are  only  tangentially  related  to  the  theme  of  Volume 
I.  Several  references  were  made  in  Volume  1  to  appendices  in  this  volume 
for  the  details  of  specific  methods. 

The  notation  introduced  in  each  appendix  applies  only  to  that  append!^ 
unless  otherwise  noted. 
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least  squares  solution  to  U.6)  will  give  biased  estimates. 
When  we  use  leant  squares  to  minimise  che  residuals  by 

rf-  X  u*2‘°  • 

m  K^O 

ve  obtain  the  final  expression 


(A. 7) 


(A.  8) 


N-l  Y-l  Y-l 

plo  “p  Jo  ?P+K  ‘  "  Jo  Yn+k 


We  will  get  biased  estimates  for  the  a  . 

P 


m*0, . . . ,N-i 


(A. 9) 


One  way  to  correct  this  Is  to  use  a  method  known  as  Iterative 
generalized  least  squares.  Rewrite  (A. 7)  as 


X  a_  m  *  K"0»  1» « « •  »Y“1  « 

pfo  P  P+K  K 


(A.  8) 


Now  define  new  notation  for  the  sake  of  convenience.  First  Introduce  the  shift 
operator  q  defined  so  that 


qf  ■  f 
q  K  K+l 


q  f  m  f 

q  K  K+2 


(A.  9) 


The  polynomial  operator  A(q)  Is  defined  as 
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A(q)  -  j  a  q® 
m-0  m 


(A. 10) 


A- 2 
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Now  (A. 8)  can  be  written  as 


A(q)  Yr  -  WR  ,  let  a  -  a  .  (A. 11) 

In  matrix  form  (A. 8)  is 


“Y„  Y,  Y„1 

~a7 

rw,i 
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and  (A.  11)  is 
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Now  let  us  operate  on  WR  with  the  polynomial  filter 

B(q)  \  m  \  (A. 12) 

to  give  nR  which  is  an  uncorrelated  noise  sequence  (random  variables)  so  that 
B(q)  A(q)  YR  -  . 

Let  A  commute  with  B  so  that 


A(q)  B(q)  YR  -  nR 


(A. 13) 


A- 3 


and  define 


\  -  B(q)  Yx  (A. 14) 

so  that 

A(q)  Yr  -  nR  .  (A. 15) 

Now  equation  (A. 15)  has  uncorrelated  residuals  and  hence  the  solution  should 
converge  to  unbiased  estimates. 

The  iterative  procedure  is  therefore  as  follows: 

1.  Solve  equation  (A. 11)  using  the  normal  least  squares  procedure. 

This  gives  an  estimate  of  A(q)  which  can  be  thought  of  as  the 

A 

ci's  of  (A. 8).  Call  that  estimate  A(q) . 

A 

2.  Substitute  A(q)  into  (A. 11)  to  produce  an  estimate  of  the 

A 

residuals  WR. 

A  A 

3.  Use  Wg  in  (A. 12)  to  make  a  least  squares  estimate  of  B(q) . 

A 

4.  Calculate  Yr  in  (A. 14). 

A 

5.  Use  the  YR  to  make  a  new  least  squares  estimate 

6.  Continue  this  procedure  until 

T— *  2 

2,  W  no  longer  decreases  with  the  next  iteration. 

p-0  K 

The  question  now  is: 

A 

How  does  one  obtain  B(q)  of  step  3  above? 

Assume  we  have  used  the  least  squares  process  once  to  find  the  coefficients  a  . 

P 


A-4 


Next  calculate  the  y  residuals  W  as 

& 


N 

Z  ao  ^rvfK  "  *  K-0,1, . . .  ,y-l 

p-o  p  P+K  K 


The  results  of  steps  1  and  2  above  gives  a  y  diaenHional  vector  of  the  re¬ 
siduals  w  . 

Calculate  the  autovariance  function  of  the  by 
r(u)  «  ,  u-0,1 . y-1 


where 


1  Y-u 

c(u)  -  *  I  (W1  -  W)  <W1+u  -  W)  ,  u-0,1 . y-1 


-  1  i 

W  -  *  t  W  . 

Y  i-1  1 


Hence  the  autocovariance  function  r(u)  is  defined 

Y  («t  -  5;  <wt+u  -  5) 

r(u)  -  - -  ,  y-0,1, . . .  ,y-l 

V  <w4  -  W) 2 

i-i  1 


We  now  have  a  vector  R,  y  long. 


We  want  to  test  this  R  vector  for  whiteness  of  the  residuals  W^, 
The  standard  deviation  of  a  single  autocovariance  function  estimate  is 


. . -  iw.i  Mr' '  r«  “1 1  U  ~  1’i  t'lVl  f'  *  *  •'* 


The  95%  confidence  limits  for  a  single  autoccvariance  sample  are 


r(u)  +  i,96 

ST 

1  96 

Now  supposedly  if  95%  of  the  samples  r(u)  are  less  than  ■=  —  then 

Sy 

we  have  95%  confidence  that  the  noise  is  white. 

If  the  r(u)  flunks  the  whiteness  test  which  will  happen  the  first  few  times 
we  Iterate  we  must  then  go  back  and  whiten  the  residuals  WL. 

From  Equation  (A. 12)  we  want 

B(,)  wr  *  \  • 

L 

Remembering  that  B(q)  -  J  b  qm  , 

m«0 


As  an  example,  assume  L  ■  2  and  K  *  3.  Then  we  get 


W0  wi  w2 

W1  W2  W3 

w2  W3  W4 

Wo  W.  W, 


Solving  Wb  -  n  for  the  bQ  by  using  least  squares  gives 


T  T 

WAWb  -  W  n 


b  -  (WTW)  WYn  . 


Therefore  the  new  Y_  previously  defined  as 


can  be  written  as 


h  ■  Vk  +  Vk+i  +  V*«  '  +  V: 


K+L 


APPENDIX  B 

EFFECT  OF  IN^KEASING  MODEL  ORDER  IN  PRONY'S  METHOD 


A  study  of  the  effect  of  increasing  model  order  in  Prony's  method 
has  been  made  and  the  results  are  presented  here  along  with  aome  preliminary 
conclusions.  In  past  investigations  of  Prony's  method  it  has  beeu  noted 
that  increasing  the  model  order  above  the  known  order  of  the  waveform 
improves  Prony's  method's  ability  to  estimate  accurately  the  true  poles. 
Although  the  accuracy  of  the  pole  estimates  is  improved,  a  side  effect 
of  this  procedure  is  the  problem  of  distinguishing  between  true  poles 
and  those  poles  that  simply  fit  to  the  noise  and  have  no  relation  to  the 
Information  that  is  to  be  extracted  from  the  waveform. 

In  this  study  we  relate  the  inaccuracy  or  bias  of  the  parameter 
estimates  to  the  amplitude  of  the  residuals  of  the  least-square  Prony 
procedure.  We  assume  that  the  inhomogenous  solution  (defined  in  Volume 
I,  Section  2)  is  used  to  find  a  parameter  vector.  The  term  "residuals" 
is  Identical  to  "equation  error".  In  this  appendix,  N  is  the  number  of 
poles  modeled  and  M  is  the  number  of  samples.  When  the  residuals  are 
zero  the  least-squares  Prony's  method  either  has  reduced  to  curve-fitting 
Prony's  method  (M-2N)  or  is  processing  a  noise-free  waveform  at  the  proper 
model  order.  In  this  case,  the  pole  estimates  are  quite  accurate  and  are 
free  of  bias.  Conversely  the  bias  in  Prony’s  method  is  directly  related 
to  the  magnitude  of  i  residuals. 

Tables  B.l  through  B.9  display  the  effect  of  increasing  model  order 

at  three  different  noise  levels.  is  the  standard  deviation  of  the 

N 

noise,  is  the  average  standard  deviation  of  the  residuals  over  ten 
Monte  Carlo  runs.  Each  table  shows  the  true  parameter  values  in  the  first 
column,  the  average  parameter  values  over  ten  Mone  Carlo  runs,  and  the 
variance  of  each  parameter  in  the  third  column.  For  all  nine  cases  the 
time  window  size  is  kept  constant.  Thus,  M  ■  100  and  AT  ■  0.13  seconds. 


B-l 


The  chree  pole  pairs  shown  under  the  average  coition  were  obtained  by 
searching  all  N  poles  to  find  the  pairs  that  were  closest  to  the  true 
poles  In  value.  In  general,  of  the  extra  poles  not  reported  here,  N-6, 
have  residues  two  orders  of  magnitude  below  those  of  the  true  poles. 

From  the  results  of  Tables  B,1  through  B.9,  the  following  observations 
can  be  made: 

1.  Prony’s  method  seems  to  guarantee  good  results  If  oR  <  a^. 

2.  an  <  a.,  seems  to  occur  when  4N  >  M. 

That  good  results  begin  at  oR  ■  c?^  Is  not  surprising.  In  this  case  one 
would  expect  that  the  residuals  are  beginning  to  approximate  the  noise  and 
that  the  "modeled"  waveform  approaches  the  uncorrupted  waveform.  But  the 
observation  that  cR  ■  approximately  when  4n  ■  M  has  no  obvious  explan¬ 
ation.  This  observation  would  obviously  not  hold  true  if  we  were  to  apply 
Prony’s  method  to  a  noise-free  waveform.  In  this  case,  the  residuals 
would  drop  to  zero  as  soon  as  the  model  order  reached  or  exceeded  the  order 
of  the  waveform.  We  might  surmise  than  that  thir  occurranca  at  4N  -  M 
is  attributable  to  the  nature  of  the  uncorrelated  noise. 

There  is  no  simple  explanation  for  improved  accuracy  at  higher 
orders.  At  least  two  factors  seem  to  be  involved: 

1.  Increasing  the  length  of  the  parameter  vector  (by  Increasing 
the  order)  is  an  effective  treatment  for  the  dense  sampling 
problem  (which  is  described  in  Volume  I,  Section  3  of  this 
report.  ) 

2.  Making  the  data  matrix  more  nearly  square  must  reduce  the 
magnitude  of  the  residuals.  Smaller  residuals  Implies  smaller 
bias  in  the  parameters. 
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Table  3.1. 


Model  order  study  M-100,  N-12,  o  -. 
5R-.18O0.  N 


100, 


TRUE  AVERAGE 


-8.200000E-02  ‘-6.  8  129081-01 - 7* 946325r-03' 

> 8*200000 E- 02  -6.  812908E-01  7.546325E-03 
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-1.470000E-01  '  -6*  64202  8E-Q1 - 3.8 21274T*03‘ 

•1* 880000E-01  -3*  329218E-01  2.371370E-03 

•1. 88 0 00 0E-01  -3*  329216E-0 1  2.371370E-03 


“9726  0  0  00  E-01  0*  1 - 07 - 

•9. 26 OOQOE-Ol  0*  0* 

2. 8740Q0E+0Q  2.  341076E«-00  3.239736E-03 

- 2. 8740 00 E ♦  00 -2.  3410768 ♦  00 - STHWCTCTJT 

4.839000E400  4. 84431 1E*00  6.047468E-04 

-4.835000E+00  *4.  8443UE»00  6.047468E-04 


'  l.OOOCOOEfOO . 1.  1093I9E«-00 - 779'0'9roE=TT2' 

1*  QQOOOOE+OO  1* 1Q9399E+00  7.509193E-02 

l.QQOOOOE+OO  1*14437  8E  +  00  1* 381987E-02 

*Ti  0  0  0  00  0  E*  0  0 - 1.T4437  8EVOO - lT3’mT7E=^^^• 

1*  QOQOOOE4  00  1*  189903E+00  6.474729E-03 

1*  OOQOOOE+QO  1,  189903E»00  6.474729E-03 


- 4*4W111E-“15 - 47077TOE»2B- 

0*  4*  4 19111E-15  4*  0  776  92E-28 

0*  1.907967E-01  1.508953E-02 

0* '  -1*  5  07967E-01 f75W«E^02 


0*  -1*  701996E-01  9* 4768 87E-03 

0 .  . .  1.  701996E-01  9*4768  87E-03 


REAL  PARTS 
OF  POLES 


IMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 
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Table  B.2. 


Model  order  study  M«100»  N«20,  oN«.100v 
aR".l24. 


TRUE 

'T.200000E-U2 
•8.200000E-02 
■l*4700Q0E»Ql 
> 1*470 000 £*01 
■1.88 OOOOE-Ol 
■1.880000E-01 


AVERAGE 

■9. 069350E-02 
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4.83900QE+00  4.8402S1E*00 

-4.833000E+00  -4. 840 21 1E*00 


3.637123E-05 
3.657123E-03 
1*  011044E-0V 

TTiTfnOT=Tf 

1.176039E-04 

1.176039C-04 


IMAG.  PARTS 
OF  POLES 


l.  0 0  0  0 do E ♦' 0 0  '  “ ' 7  U  0‘27 8 5 9 E ♦  0 0  2.  0  3744 3 E-  03 

1*  OOOOOOE  +  OO  1*  0  278S9E  +  00  2.037443E-03 

„l*m0Q0  E  too. _ l*0.?Afk&71±Jlll _ 2iliaiA9  fcfli 

1. 0  0QQQ0E  +  00  1*  0364S7E  +  00  2.199189E-03 

l.OOQOOOE+QO  1*  0 18 026E  +  00  1.059788E-03 

1. QOQOOOE+QO _ u  0 18  028E*  DO  1.  0 397_88E-03 


MAGNITUDES 
OF  RESIDUES 


O'* 

-1.230493E-02 

9.62Q293E-04 

0* 

1. 230493E-02 

9.820293E-04 

0* 

-1.928778E-02 

2.042011E-03 

0* 

l.^SttOE^OE 

2.0420 11E-03 

0* 

-E.326322E-03 

2.  Q243  03E-Q3 

0* 

5.328322E-03 

2.024303E-03 

RADIAN  PHASE 
OF  RESIDUES 
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Table  B.3.  Model  order  study  M-100,  N-36,  o„-.100, 
3„-.0678. 


TRUE 

8.200000E-02 
6*  2000002*02 
1.470000E-01 
1.470000E-01 
1.8800002-01 
1.8800002-01 


AVERAGE  Q 

-6*.  16917 5*2-02  3.  0566572-05 

•6. 1691712-02  3.0566572-03 

. .  -1. >5353.52-01 _ VO7.227JE.-0> 

-1.4835382-01  1.1722732-04 

-1.0425902-01  1.8613972-04 

._=0fJ»mUbll _ UM1W1-04 


REAL  PARTS 
OF  POLES 


9.  26 0  0 OOf- Oi" 
9. 260000E-01 
2. 6740002+00 
2. 8740002+00 
4. 6350002+00 
4. 6350002+00 


•9.2664802-01  2.3949362-03 

2. 8748362+00  5.3018692-05 

•  2.8  n  a  m  ♦  on - 573*1)1 8OT2-  0  i 

4. 8361282+00  4.4500362-05 

•4.83612SC+0Q  4.4500362-03 


IMAG.  PARTS 
07  POLES 


1. 0000002*00 
1. 0000002*00 
1. 0000002+09 
1.  0000002+00 
1. 0000002+00 
1. 0000002+08 


9.9310762-01  1.3250782-03 
9.9310782-01  1.3250782-03 
.1. 0137202+00  4.1196332-03 
1.0137202+00  4.1196332-03 
9.6910632-01  6.2767592-03 
9.  6 91 0 632-01 _ 6.2767592-03 


MAGNITUDES 
07  RESIDUES 


•3.  *563 1472-03 
3.5631472-03 
1.  4453582-03 
•  1.445*582=11 
•1.4026732-02 
1.  4026732-02 


8.5535592-04 

8.5535592-04 

IjilffttibJUL 

177902252-03 

7.1345402-04 

7.1345402-04 


RADIAN  PHASE 
07  RESIDUES 


Tfebla  B.4.  Modal  ordar  itudy  M-100,  N-8,  c^-. 00100, 
oR-. 007098. 


TRUE 

a  NBA.  flUA  -A  ~T  M  AM  »  IBM 

AVERAGE 

a 

tnim -TT* 
-8.2000008-02 
-1. 4700006-01 

-5774682)8-01 
-4,  7460212-01 
-3.23707)8-01 

^.f7l5TO?33 
6.  6719  06E-03 
1. 1543838-03 

•1.  47Q0fltJE-()i 

-l.seooooE-oi 

-1.880000E-01 

-3.217'07nr“ffT" 

-1.9949998-01 

-1.9949191-01 

“Tm'4irnnr 

2. 614657E-P5 
2. 614697E-0) 

• 

9.2600008-01 
-9.260000E-01 
2* 874000T4  0-0 . 

8.99030)8-01 

-6.99030)8-01 

2. 8743948+00 

8. 7045 12E-04 
6. 704912E-04 
2.6016898-04 

-i.iHoooeioo 
4*  839 0008400 
-4.8350008400 

-2. 874394I+00 
4.6432412400 
•4.8432418400 

2.  6016892-04 

3. f 698662-0) 
3. 669866E-0) 

1«  0 u uOOu  c  □ 
1.000000E+QO 
1.0000002400 

177371898400 

1.7371998400 

1.2238202400 

9.1143 46E-0T 
9.1143668-03 
1.4617208-03 

l*  ooooooEioil 
It  OOOOOOE  +  OQ 
1. OQOOOOE+QO 

raimrur 

9. 0263708-01 

9.  02837  0e-01 

1. 4)l720K-03 
4. 5900  50E-04 
4.950050E-04 

REAL  PARTS 
OF  POLES 


IMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


"0.- 

2.  0  228  031-01"“ 

“77ZT5594T-03 

0. 

-2.0228038-01 

7.2359948-03 

0. 

2.29633)2-01 

3. 769424E-03 

0. 

-2. 2983j)X-01 

2.7694248-03 

0. 

1.2808948-01 

2.4557898-05 

0. 

-1.2606948-01 

2.49576IE-0S 

RADIAN  PHASE 
OF  RESIDUES 


Table  B.5.  Model  order  study  M-1G0,  N-12,  a 


TRUE 

3r-. 00300. 

AVERAGE 

A 

-8. 200000E-02 

-6.3347606-02 

3. S19636E-Q5 

•6. 200000E-02 

•6. 3347S0E-Q  2 

3.819636E-06 

-.1  *.47.0  0  001-0 1 _ - 1 1.482  0  4i  E-.0i_. 

_3*. 3.0  75561-06. 

-1. 470000E-01 

-1.  482045  E-01 

3. 30  75  56E-QS 

•1. 66  OOOOE-Ol 

-t.  679567  E-01 

1. 655490E-07 

| 

I 

■ 

,.JkjJUU&9UL6±0JL 

9.260000E-01  9. 2632276-01 

-9.260000E-01  -9. 2632276-01 

if*  835000E+00  4.  6390411+00 

-4*  6390Q0E+Q0  -4. S3904II+00 


3. 0542061-06 
3. 0  542  06E-06 
2.313336E-06 

iltimaHi 

6.7570916-07 
6*  7570516-07 


1. OOOOOOE+OO  1. 00544L6+00  9.2241066-05 

1. 00  Q0QQE  +  0Q  1.0054416+00  9.2241066-05 

^l.aMflJaEtPQ  .,  4.7-160  98 W1 

1. OOOOOOE+OO  1.0038536+00  4. 7160966-05 

1. OOOOOOE+OO  9. 994317E-01  4.064030E-06 

l.OQOOQflg+QO  .9.  994  3171-01  4.  0640  306-0  Bl 


0* 

-1.1104116-03 

6. 015323E-05 

0. 

1. 11041} 6-03 

8. 0153236-05 

0. 

6.  180  273 1-04 

4.066317E-0S 

0. 

-8.1802706-04 

4. 066317E-05 

0. 

1.6991516-03 

9. 11 37 11E-06 

_1* _ 

-1.  lilUIMl 

...  9.1137111-06, 
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-.00100, 


REAL  PARTS 
OF  POLES 


IMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 


Table  B.6.  Model  order  study  M-1G0,  N-20,  aN-.001, 
crR-.  001246. 


TRUE 

•8.200000E-02 
*8*20  00006-02 


AVERAGE 

'8.  199893E-Q2 
•8.  199895E-02 


6.9594356-09 

6.939433E-09 


wuw w  * 

-1.470000E-01 
-1.88OOO0E-O1 
■li 98  0QQ0E-Q1 

V.  Zl  *****  y  * 

-1.47015SE-01 

-1.8303S8E-01 

-1.880399E-01 

1.0346476-05 

1.1212236-08 

9.2600Q0E-01 
-9. 260000E-01 
2.8740QQ£*00 

9.2600796-01 

-9.2600798-01 

2*  B7402IE*00 

3.504264E-03 

3.304264E-03 

9.5360526-09 

-2.8740006*00 
4*  835000E* 00 
-4.-833Q_0P_L»M_ 

-2.  874025H-00 

4.  8349946*00 
-4.83499*6*00 

9.536032E-09 
1.374833E-08 
1«  3745336-09 

1.0000006*03 
1.0000006*00 
1*  0000006*00 

9.9995096-01 

9*  9995396-01 
1.0001106*00 

2* 1308856-07 
2.1308536-07 
_  2.J.91616E-Q7 

i.OOQOOOE*OQ 
1*  QOOOOOE>QO 
1.  0000006*00 

1*  0001106*00 

1. 0000666*00 

1.  00006Se*00 

2.0916146-07 
9*  0947372*08 
9.0947376-08 

0. 

0. 

-4*  4 011026-05 

4*  40190 22*09 
-8. 9732606-05 

9. Q35431E-T8 
9.0354316-08 
2.  Q  045726-07 

0. 

0. 

JU _ 

6. 9732906-05 
-2. 9039196-03 
2.9Q39UE-05- 

2.  0049726-07 
2. 2595966-07 
-.Jiili ZSS&zAL 

REAL  FARTS 
OF  POLES 


XMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 
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Table  B.7.  Model  order  study  M-100,  N-8,  O..-.0100, 
or-.0394. 


TRUE  AVERAGE  a2 

8/20  00 00 £-02  -4,  8 1594  IE-01  4.313S74E-Q3 

8«  20  Q000E-Q2  -4. 815941E-01  4.313874E-03 

1.470000E-01  -4« 931 019  E-Ql  1.633214E-03 

1.470000E-01  -4.931013E-Q1  1.6332UE-03 

1.8800005-01  -2*  503945E-Q1  4.518706E-05 


REAL  PARTS 
OF  POLES 


9.2600005-01 
9.260000E-01 
2.87400Qg»0Q 
2.874000E4-00 
4* 83S000E+  00 


■2.44487SE+0O 

4.820967E+00 


6.292900E-04 
1.6992 93E-04 


IMAG.  PARTS 
OF  POLES 


l.OOOOQOS^OO  8.254609E-01  2.930682E-Q2 
1*  Q00080E+  00  8* 254609E-01  2.930682E-02 
1. 00  OOOQE»OQ  1* 154Q58E+Q0  3.968307E-03 
1. 00 0000 **00  l.tiVoV8E400  3.S58307E-Q3 
l.OOOOOOEfOO  li 1 12501 E^00  4.643365E-Q5 
1» 0000005*00  1. 11230 1E+00  4.643363E-05 


MAGNITUDES 
OF  RESIDUES 


2*  991933E-15 
1»  623  083  E-01 
•1.6230S3E-01 
•1.  075143  E-01 


6*  5  963*»lE-28 
2.  7298  86E-03 
2.729888E-03 
9.  082337E-04 


RADIAN  PHASE 
OF  RESIDUES 


'l 


I 


I 


ia 


T^ble  B.8.  Model  order  study  M-100,  N-12,  g  «.0100, 
3  0295. 


TRUE  AVERAGE  g2 

8 . 2  0  0  0  0  0  E -02  -2 . 35316iE-fll  2a  <3592  68E-03 

8.  2000005-02  -2. 35315  IE- 0 1  2.859268E-03 
1. 470 POPE- 01  -2.  4074SSE-Q1  5.513517E-QV 

1. 470000E-01  -2. 407455E-Q1  5. 5 1 8517E-0'. 

l.anOOOOE-Ql  -1.9345&1S-01  1.094043E-05 

1 . SHO000E-Q1  -1. 934551E-01  1.094Q43E-Q5 


REAL  PARTS 
OF  POLES 


9. 26  0  00  0S- 01 
2.874000E+Q0 

9.4d447J£-di 
-9. 408873E-01 
2.  899695  E+  00 

JTT5TT53^ir 

3.593158E-04 

1.687968E-04 

^2737^  000  ETO  fl~ 
4.835000E+00 
•4,8390005+00 

rnnr 

4, 842363  E+00 
-4, 842353E+00 

r.",6'87g,6,4i-u; 
1, 1226 18E-Q4 
1.122S18E-Q4 

1.  Q0OQ0QE+QQ 
1.000000E+00 
1*  OOOOOOE+OQ 

1. 4i60SQE+00 
1*  426 05 0 E+0  0 
1.211075 E+00 

1.70  80  38E-02 
v  1.70 80 3 BE- 02 
v*  2, 967772E-03 

1,  OOOOOOE+OQ 
1. OOOOOOE+QO 
1. OOOOOOE+OQ 

1. 211 075E+00 
9. 755319E-01 
9.795319E-01 

2.967772E-03 
4, 1906  53E-04 
4.190653E-04 

s.  nmye-ai — Tttzrunn 

2.  7 12845E-03  7.626122E-03 
7*  308093E-Q2  3.766332E-03 

77TJnTT5T^ff? - 3;76T332E-d3 

1.0259S3E-Q1  8.128753E-0V 

1.  025580E-01  8*1287  53E-04 


IMAG.  PARTS 
OF  POLES 


MAGNITUDE 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 


1 


i 

I 
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APPENDIX  C 

THE  ADAPTIVE  METHOD  FOR  RESONANCE ESTIMATION 

The  adaptive  method  reeulta  from  the  generalized  scheme  of  Volume  I* 
Section  2,  under  the  following  assumptions : 

1.  a  “  1  and  g  *0. 

o  o 

2.  Pj  ■  — - —  ,  i  ■  l,...n  and  P  ■  1. 

i  z-z^  *  o 

3.  The  model  input  is  a  unit  sample  at  k  -  0  (discrete  impulse). 


i 

} 


The  unique  feature  of  the  method  is  that  the  filter  poles,  z^,  may  be 
adjusted  to  any  value  in  the  Z-plane,  An  adaptive  technique  fnr  adjusting 
the  filters  consists  of  first  initializing  the  filters  to  arbitrary  values 
in  the  Z-plane  and  repeating  the  following  steps: 

1.  Find  an  estimate  of  the  process  transfer  function  using 
the  current  filters  in  the  model. 

2.  Set  each  filter  pole  to  one  pole  of  the  estimated 
transfer  function. 

This  procedure  is  repeated  until  a1  approach  zero.  The  poles  of  the 
process  can  be  estimated  during  the  course  of  iteration  by 

A  *i 

zi  1  +  at 


removing  the  need  for  finding  the  roots  of  a  polynomial.  The  filter  poles 

are  updated  to  on  each  iteration.  When  the  a^  approach  zero  the  pole 

updating  ceases  and  the  method  converges.  At  each  iteration,  the  a.  are 

1 

found  using  any  of  the  techniques  for  finding  a  parameter  vector  described 
in  Volume  I,  Section  2.  Perhaps  the  simplest  method  is  to  choose  the 
parameter  vector  as  the  weakest  eigenvector  n f  ft*Q.  At  convert  jnce  the 
S-plane  poles,  can  be  obtained  from  the  fi’.tsr  poles  by 


s 


i 


in  z^ 
T 


and  the  A^  •  3^ 
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The  adaptive  filtering  method  seems  to  be  a  very  useful  technique 
because  it  uses  what  seems  to  be  optimal  filters,  band  pass  filters,  for 
estimating  pole  locations.  In  addition  it  presents  a  measure  of  the  error 
and  iteratively  adapts  the  filters  until  the  error  is  at  the  level  where 
it  is  desired. 

As  an  example  of  the  use  of  the  adaptive  filtering  method  consider 
the  data  shown  in  Figure  C.l.  These  data  were  numerically  generated  by 
using  the  time  domain  computer  code  TWTD  [C.l],  The  structure  modeled 
by  TwTD  was  a  thin  cylindrical  scatterer.  The  signal  was  contaminated 
with  noise  to  give  a  15  dB  signal-to-noise  ratio.  Figure  C.2  shows  the 
resulting  poles  from  five  Monte  Carlo  trials. 

The  following  statements  can  be  made  about  the  adaptive  method  after 
studying  it. 

The  adaptive  method  is  a  new  method  which,  in  many  cases,  pro'  ides 
excellent  pole  estimates  under  difficult  conditions.  The  method  is 
unique  in  that  a  solution  to  a  polynomial  is  not  required  to  find  estimates 
of  the  process  p.les.  the  method,  in  effect,  "swallows"  the  polynomial 
solver  in  its  own  iterative  pole-searching  scheme. 

Attempts  to  analyze  waveforms  consisting  of  highly  damped  exponential 
components,  such  as  the  transient  responses  of  a  sphere,  have  not  been 
successful.  The  adaptive  method  does  not  converge  for  waveforms  which 
display  double  pole  characteristics,  that  is,  waveforms  with  components 
of  the  form  t  exp(st).  Slight  modifications  to  the  adaptive  method  might 
allow  the  analysis  of  such  waveforms. 
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APPENDIX  D 


THE  PENCIL-OF-FUNCTIONS  METHOD 

Tha  pencil-of-functions  msthod  rssults  from  ths  generalized  modsl 
describsd  in  Volums  I,  Ssction  2  with  ths  following  assumptions: 

1.  Oq  ■  1  and  8q  ■  0. 

2.  ■  F^Cs)  ■  (l/s)1  (cascadad  continuous  intagrators) . 

This  mathod  is  not  vary  assy  to  implamant  on  a  digital  computar  sinca  tha 
continuous- tirna  intagrators  cannot  ba  imp lamented  exactly  by  any  algorithm. 
Whan  approximate  integrators  ara  cascaded,  as  they  are  in  tha  pancil-of- 
f unctions  mathod,  largo  errors  can  ba  quickly  accomulated  and  tha  intended 
result  after  a  number  of  integrations  destroyed. 

This  difficulty  can  ba  resolved  by  cascading  discrete  intagrators  to 
obtain  filters  with  pulsa  transfer  functions  given  by: 

rt  *  'i(,>  ■  (A)  ■  (if 

which  can  ba  implemented  on  a  digital  computer  with  no  error  by  using 
difference  equations.  The  variable  Z  is  defined  as 


2  "  1  * 


and  z  is  the  z-transform  variable. 

Whan  discrete  integrators  are  used  the  discrete  pencil-of-functions 
mathod  results.  The  poles  of  the  pulsa  transfer  function  of  the  process 
or  waveform  may  be  estimated  as 


si  “  l-V 


kL 

where  is  ths  i  zsro  of 

*n  +  a.  Zn-1  +  ...  +  a  -  Q  (D.l) 

l  n 

In  ths  original  pencil-of-functions  mtthod  ths  wars  found  by  using 
Jsln's  msthod  for  constructing  ths  parameter  vactor  which  is  dascrlbad  in 
Voluma  I,  Saction  2.  With  Jain's  mathod 


aU  ak 

whara  A..,  in  this  csss,  is  tha  olsmant  at  tha  i  row  and  j  column  of 
»dj  fl*n.  Othar  methods  can  ba  usad  to  astimata  tha  paramatar  vactor. 

Tha  S-plane  astimataa  of  tha  poles,  s^,  ara  ralatad  to  tha  zaroa  of 
(D.l)  by 

i 

Ona  difficulty  with  tha  pencil-of-f unctions  mathod  is  ralatad  to 
tha  attenuation  of  tha  higher  frequency  modes  of  tha  process  output  by 
tha  rapaatad  integrations  applied  to  the  output  waveform.  It  can  ba 
verified  that  an  integrator  is  simply  a  first  order  filter  whose  Laplace 
transfer  function  has  a  pole  at  tha  origin  in  tha  S-plana.  Such  a  filter 
tends  to  suppress  tha  higher  frequencies  present  at  its  input.  Tha 
higher  frequency  suppression  phenomenon  is  illustrated  in  Figure  D.l. 
Normally,  when  an  exponential  function  is  integrated  repeatedly, 
components  of  power  of  time  exist  in  the  higher  integrals  as  wall  as  the 
original  exponential  function  components.  In  Figure  D.4  the  components  of 
powers  of  time  have  bean  subtracted  from  tha  integrated  waveforms  in 
order  to  make  the  attanuation  of  the  higher  modes  more  evident.  Tha 
first  waveform  is  a  hypothetical  waveform  provided  for  analysis.  Tha 
waveforms  that  follow  are  the  integrals  of  Increasing  order  of  tha  first 
waveform  and  display  the  increasing  dominance  of  the  fundamental  mode  or 
mode  of  lowest  frequency.  Further  integrations  yield  nearly  identical 
waveforms.  The  integrated  waveforms  tend  to  become  linearly  dependent 
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vava forma .  Tha  intagratad  wavaforma  tand  to  bacoma  linaarly  dapandant  at 
highar  nodal  ordara.  Tha  matrix  than  tanda  to  alngularity  and  tha 
mathod  bacomaa  unatabla.  Tha  auppraaalon  phanomanon  occura  in  both  tha 
dlacrata  and  continuous  mathoda.  In  fact,  avan  for  vary  modaat  modal 
ordara,  tha  mathod  can  bacoma  numarlcally  ill-conditionad  to  a  dagraa  that 
apacial  cara  muat  ba  takan  to  aaaura  accurata  invarslon  of  n*n. 
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METHOD  OF  REDUNDANT  AVERAGING 


Tha  redundant  avaraglni  schema  la  a  preprocessing  schema  Chat  attampta 
Co  oomblna  redundant  data  (that  la,  more  data  than  ara  necasaary  to  determine 
tha  parana tara)  In  a  way  that  avolda  tha  blaa  lntroducad  by  ualng  a  laaat- 
aquaraa  achama  to  comblna  radundant  data.  Tha  mathod  attampta  to  tranaform 
a  raw  wavaform  of  aora  than  2N  aamplaa  whara  N  la  tha  modal  ordar  Into  a 
praprocaaaad  wavaform  of  exactly  2N  aamplaa  by  avaraglng  within  tha  wavaform. 
Tha  avaraglng  can  ba  dona  ao  that  tha  axpactatlona  of  tha  polaa  of  tha  pra¬ 
procaaaad  wavaform  ara  aqulvalant  to  tha  axpactatlona  of  tha  polaa  of  tha 
raw  wavaform  provldad  tha  addltlva  nolaa  on  tha  raw  wavaform  la  aaro  maan 
and  uncorralatad  batwaan  sucasalva  aamplaa.  Tha  praprocaaaad  wavaform  may 
than  ba  procaaaad  with  cuxva-flttlng  Prony'a  mathod  to  avoid  tha  blaa  of  tha 
laaat-aquaraa  procadura. 

Tha  moat  ganaral  daacrlptlon  of  tha  radundant  avaraglng  achama  can 
ba  atatad  almply  aa 

v1 

'  *  jlo  Vu<,  ’ k"0’1 . 21,-1 

aU 

whara  N  la  tha  modal  ordar,  X.  danotaa  tha  k  aampla  of  tha  praprocaaaad 

wavaform,  and  y^  danotaa  tha  k  aampla  of  tha  raw  wavaform.  In  ordar  to 

limit  thla  daacrlptlon  to  tha  aaaantlal  faaturaa  of  tha  radundant  avaraglng 

achama,  N. ,  tha  numbar  of  avaragaa,  an''  N  ,  tha  dacimatlon  apoch,  ara  not 
A  ■ 

axpllcltly  dafinad  hara.  Thaaa  parameters  ara  chooaan  aa  daalrad  but  uaually 
in  a  way  that  would  produca  a  daalrabla  praprocaaalng  moda  in  soma  sanaa. 

For  inatanca,  tha  value  for  tha  dacimatlon  epoch  might  ba  choaon  on  tha 
basis  of  tha  maximum  frequency,  <u  ,  praaant  in  tha  raw  wavaform  (if  thla 
information  is  available)  and  tha  numbar  of  avaragaa  might  be  choaan  ao 
that  every  aampla  in  tha  raw  wavaform  la  used  once  tha  value  of  tha 
dacimatlon  apoch  la  sat.  Hanes  tha  sampling  interval,  At,  for  tha 
praprocaaaad  wavaform  could  ba  obtained  aa 


irfi’iir*  'iiiuaa^t 


At  - 

where  ■ 

The  number 


N  At 

•  raw 


Integer  Part 


of  avaragaa,  N^,  can  ba  computed  by 


NA  -  M  -  N-  (2N-1) 


whara  H  la  tha  number  of  aamplea  In  the  raw  waveform  and  N  la  the  order 
daalred.  It  ahould  ba  noted  that  tha  above  method  for  determining  and 
Nj  la  not  tha  only  poaalbla  technique. 


A 


Tha  redundant-averaging  procedure  la  affectively  two  operational 

1.  Apply  a  low-paaa  moving-average  filter,  A(z),  to  tha  raw  waveform. 

2.  Decimate  tha  filtered  waveform  with  decimation  epoch  Ng. 


Tha  tranafar  function  of  tha  moving-average  filter  iat 

N.-l  N.-2  K  i 

A(a)  -  a  *  +  a  A  +  ...  +  a+1  • 


The  numarator  of  A(r)  can  be  factored  into 


II  (z-z.), 
i-1  1 


The  zeroa,  z^  of  the  numerator  of  A(z)  have  unit  magnlcudaa  and  argumanta, 


erg  zt 


for  i  ■  1,  N^.  It  la  clear  that  tha  firat  factor  cancala  with  tha 

denominator  giving 
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A(z)  ■  H  (z-z.), 
i-2  1 
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Table  E.l  True  poles  and  results  of  five  Monte  Carlo  runs  for  the 
redundant  averaging  example  shown  in  Figure  E.l. 


REAL  PARTS  OF  POLES 


TRUE 

POLES 


MONTE  CARLO  RUNS 


2  ' 

3 

-.127 

-.122 

-.157 

-.239 

-.103 

-.227 

AVER. 


-.107 

-.187 

-.271 


Z  DEV. 


IMAGINARY  PARTS  OF  POLES 


MONTE  CARLO  RUNS 


TRUE 

POLES 


.926 

2.874 

4.835 


2 


.904 

2.913 

4.899 


•  .942 
2.889 
5.065 


5 


.936 

2.882 

4.803 


AVER. 


.925 

2.902 

4.899 


X  DEV. 


A(z)  chen  has  zeros  which  are  evenly  spaced  around  the  unit  circle  but 
the  zero  at  z»l  is  canceled.  Since  there  is  never  a  zero  at  z»l,  the 
moving-average  filter  can  be  viewed  as  a  type  of  low-pass  filter,  A 
z-plane  plot  of  the  zeros  is  shown  in  Figure  E.l. 

As  an  example  of  the  application  of  this  method  consider  the  noise- 
free  waveform  shown  in  Figure  E.2a.  This  waveform  consists  of  400  data 
points  representing  a  sixth  order  exponential  function  generated  with  the 
poles  listed  in  Table  E.l.  This  data  was  then  corrupted  by  adding  white 
noise  with  a  standard  deviation,  o,  of  0.5.  Figure  E.2b  shows  the  noise 
contaminated  signal.  This  signal  has  a  signal  to  noise  ratio  of  15.6  dB 
where  the  signal  to  noise  ratio  is  defined  as 

S/N  -  20  log 

where  R  is  the  peak  amplitude  of  the  transient  signal,  Five  Monte 
Carlo  trials  were  run  on  this  data  using  the  redundant  averaging  mathod. 

The  model  order  was  selected  to  be  8  (two  more  than  the  known  order)  and 
Ng  and  were  chosen  to  be  16  and  160  respectively.  Figure  3.2b  shows 
the  reconstructed  waveform  obtained  from  one  of  the  Monte  Carlo  trials. 

Note  how  good  the  fit  is  considering  the  signal  to  noise  ratio  was  15.6  dB. 
Table  E.l  lists  the  results  of  the  five  trials  and  shows  the  pole  averages, 
standard  deviations  and  per  cent  deviations.  It  should  also  be  noted  that 
the  signal  to  noise  ratio  compared  to  each  pole  residue  Is  only  6  dB. 

Hence  from  poles  with  a  6  dB  information  content  we  were  able  to  recover 
them  very  accurately. 

The  redundant  averaging  scheme  has  two  difficulties  or  limitations 
that  can  cause  the  method  to  be  lass  effective: 

1.  If  the  sampling  rate  in  the  preprocessed  waveform  is 
sufficiently  low,  the  higher  frequency  poles  can  be 
folded,  perhaps  several  times,  About  the  Nyquist 
frequency. 


E-6 


2.  The  method  can  null  modes  in  the  preprocessed  waveform 
that  exists  in  the  raw  waveform. 

The  frequency  folding  phenomenon  introduces  an  ambiguity  in  that  it 

is  not  known  how  many  times  a  particular  pole  has  been  folded  about  the 

Nyquist  frequency.  Hence,  N  possible  poles  are  introduced  for  each 

s 

extracted  pole  by  the  redundant  averaging  scheme  where 

N  » 

3 


(AT)  preprocessed 
(AT)  raw 


Of  course,  if  the  highest  frequency  mode  of  the  waveform  is  known  to  be 
lower  in  frequency  than  the  preprocessed  waveform's  Nyquist  frequency,  the 
ambiguity  is  resolved  but,  in  general,  this  will  not  be  the  case. 

The  mode  nulling  can  occur  if  the  averaging  parameters  are  such  that 
the  zeros  of  the  resulting  low  pass  fi-ter  cancel  poles  of  the  signal 
Itself.  That  is,  if  we  define  the  preprocessed  waveform  P(z)  as 

M 

P(z)  ■  A(z)  D(z) 

resulting  from  the  averaging  process  A(z)  operating  on  the  signal  D(z), 
then  if  A(z)  and  D(z)  each  have  terms  in  common  they  can  tend  to  cancel 
each  other. 
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COLUMN  PRONY'S  METHOD 

la  this  method  the  z-plane  estimates  of  the  poles  are  the  roots  of 
the  polynomial 

„  ,  N.  ,  ,  N.2  .  >  N.N  . 

o'q  +  ct^Cz  )  +  a2(z  )  +  ...  +  a^U  )  -  0 

The  are  determined  from  the  system, 


where  ct^  is  assumed  to  be  one  and  {y^,  y^,  y^2+jj_2_}  denotes  the 

sequence  for  the  waveform  containing  N^+N  samples.  It  should  also  be 
noted  that  the.  polynomial  has  N^  roots  only  N  of  which  are  related  to 
estimates  of  the  true  constituent  poles  of  the  waveform. 

2 

Column  Prony's  method  offers  n  means  of  combining  N  +N  data  points 
compared  to  the  2N  data  points  of  the  standard  Prony's  method.  The 
column  Prony's  method  can  either  be  used  in  the  least-squares  version  or 
the  curve-fitting  version.  If  the  curve-fitting  version  of  column  Prony's 

method  is  used,  the  resulting  parameter  estimates  are  unbiased. 

2 

Unfortunately  the  method  yields  N  pole  estimates  only  N  of  which  are  the 
true  poles.  It  appears  then  that  column  Prony's  method  leads  to  the  same 
problem  that  increasing  the  model  order  led  to  in  the  standard  Prony's 
method:  an  ambiguity  in  the  identity  of  the  true  poles.  Moreover,  this 
method  has  an  additional  problem:  the  matrix  can  become  nearly  singular 
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» 


even  if  the  model  order  is  lower  than  or  equal  to  true  order.  This 
phenomenon  is  similar  to  the  singularity  of  the  standard  Prony's  method 
when  the  model  order  equals  the  waveform  order  and  the  highest  frequency 
is  equal  or  nearly  equal  to  the  Nyquist  frequency. 
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EVAN'S  AND  FISCHL'S  METHOD 


I 


» 


» 


► 
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In  Evan's  and  Fischl'a  method  [G.l]  they  define  a  "true  error" 
sequence  {e^,  k*0,  . M-l}  as  the  error  between  the  given  waveform  and 
the  "fitted"  waveform.  They  then  proceed  to  define  the  "equation  error" 

sequence  {d^,  k*0 . M-l}.  They  further  proceed  to  define  a  relation 

between  the  equation  error  and  the  true  error: 

e  ■  Wd 


where 


o*i  Vi] 
i 3  tdo  di  •“  dM-N] 


and  W  ■  A[ATA]  ^ 


where 


r*  a. 


“N  aN-l 


UO 


...  0  “1 


•  •  e  • 


0 

a0 


■  •  •  I 


aN-l 
°N  - 


M  is  the  number  of  samples,  N  is  the  number  of  poles,  a^*l,  and  the  a^; 
are  the  coefficients  of  Prony's  difference  equation  defined  in  Section  2, 
Volume  I. 


By  ueing  this  relation  between  the  two  errors  they  describe  two  iterative 
procedures  aimed  at  minimizing  tha  "true"  error  criterion: 


Ic  can  be  shown  that 


W(a)  - 


3a 


N-k 


A(a) 


-  W(a) 


3a 


A(a) 


N-k 


A(a) 


+  A(a) 


3a 


N-k 


Av  .) 


A(a)4  A(a) 


and  [(3/3aN-lt)  A(a)]  is  simply  th*  matrix  A(a)  with  onaa  raplacing  tha  aN_^ 
with  all  othar  alamants  saro. 

Tha  following  observations  can  ba  mad a  about  tha  mathodt 


1.  It  produces  optimal  pole  estimates  in  tha  sense  that  it 
minimizes  "true  error".  This  means  that  tha  mean-square  error  between  tha 
given  wavaform  and  tha  fitted  waveform  is  minimized.  These  optimal  pole 
estimates  are  obtained  only  in  the  second  iteration  phase. 

T 

2.  .ue  matrix  A  A  must  be  inverted  on  each  iteration  of  both  tha 

T 

first  and  second  procedures.  When  the  w avr.  .  u  has  over  100  samples,  A  A 
becomes  very  large,  and  hence,  expansive  tu  invert.  Consequently,  the 
method  is  prohibitively  expensive  when  the  wavaform  has  over  100  samples. 
This  method  has  yet  to  be  evaluated  by  tests  on  noisy  data.  Nothing  is 
presently  known  about  its  convergence  characteristics  of  its  tolerance  of 
noise. 


The  following  example  illustrates  the  optimal  estimates. 
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Exaaola 


N  •  1,  M  -  3 

2"  [2  1  2]  (vavaforo  eo  ba  approxiaaead ) 

s  -  lju-  t"jl 

Initial  paramatarai  ,  a  ■-! ,  •  Thasa  ara  cha  optimal  paramatara, 

tharafora  Cha  mathod  should  convarga  inanadlaCaly . 
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[3/3  3/3] 
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“  (i)  (*S)  "  “1 

Therefore  the  updated  valua  of  la  equivalent  to  ita  old  value,  hence  the 
method  hee  converged  at  the  optimal  parameters  for  the  second  Iteration 

phase. 
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APPENDIX  H 

CONSTRAINED  PRONY 'a  METHOD 
H.l  Uni  of  a  Constrained  Prony  Method 

Many  timaa,  In  tha  application  of  tha  Prony  algorithm  aoma  of  tha 
polaa  ara  known  a  priori.  It  would  ba  vary  uaaful  if  tha  Prony  algorithm 
could  ba  conatralnad  in  auch  a  mannar  that  tha  knowladga  of  tha  known 
polaa  la  uaad  in  axtracting  tha  unknown  polaa.  Tha  known  polaa  could  ba 
polaa  of  tha  driving  function  which  ara  known  from  knowladga  of  tha  Laplaca 
tranaform  of  tha  driving  function,  ayatam  polaa  which  ara  known  from  pravloua 
Prony  analyala  or  othar  tachniquaa,  or  polaa  introduced  to  modal  tha  nolaa. 

It  la  wall  known  that  for  cartaln  data  aata  Prony' a  algorithm  haa 
aotua  difficulty  In  axtracting  tha  trua  polaa  that  ara  containad  in  tha 
data.  Thla  problem  ia  ganarally  ralatad  to  tha  nolaa  in  tha  data  but  will 
not  ba  dlacuaaad  hara.  Any  mathod  by  which  tha  accuracy  of  tha  trua  polaa 
can  ba  incraaaad  will  ba  vary  uaaful. 

Logic  would  tall  ona  that  making  uaa  of  known  information  ahould 
incraaaa  tha  accuracy  of  a  calculation.  Hanca,  if  uaa  la  mada  of  known 
polaa  ln  tha  data  it  would  aaam  raaaonabla  to  aaauma  that  aoma  of  tha 
lnatabillty  ln  Prony 'a  mathod  would  ba  allavlatad.  Tha  proof  of  tha 
validity  of  thla  atatamant  raata  on  tha  actual  lmplamantatlon  of  tha 
conatralnad  Prony  algorithm  and  tha  raaulta  comparad  for  tha  aama  data 
analyzad  by  tha  unconatrainad  mathod. 

In  addition  to  aiding  ln  lncraaaing  tha  accuracy  of  tha  trua  polaa 
it  might  ba  poaaibla  to  introduca  random  polaa  which  will  modal  tha  nolaa. 

In  that  mannar,  tha  polaa  that  ara  knoim  a  priori  ara  not  tha  polaa  which 
wa  aaak.  Tha  difficulty  with  thla,  if  it  ahould  work,  ia  that  wa  naad  to 
know  tha  rank  of  tha  ayatam  ao  that  wa  can  introduca  tha  propar  numbar  of 
nolaa  polaa. 


Another  asset  of  having  a  contrained  Prony's  method  Is  that  it  can  be 
used  to  aid  in  deconvolution.  The  following  discussion  follows  directly 
from  work  performed  in  Reference  [A. 2]. 


Assume  that  a  system  is  excited  by  a  driving  function  which  can  be 
represented  in  the  form 


M 

G(t)  -  2 
j-1 


8j 


u(t). 


(H.l) 


It  io  presumed  here  that  the  driving  function  poles  s^  and  the  associated 
residues  g^  are  known.  These  can  be  determined  analytically  if  the  ana¬ 
lytical  form  of  the  driver  is  known  or  can  be  determined  from  a  Prony's 
method  fit  to  measurement  data  of  the  driving  function. 


Now  assume  that  the  response  of  the  system  to  the  driving  function 
G(t)  is 

M+N  st 

R(t)  -  T  r  e  1  u(t)  (H.2) 

i-1  1 

and  that  the  impulse  response  of  the  system  is 
N  s  t 

H(e)  -  I  h  e  k  .  (H.3) 

k-1  * 

Hence,  there  are  N  system  poles  and  residues  and  M  poles  in  the  driving 
function.  Of  course,  the  response  function  R(t)  can  be  written  as  the  con¬ 
volution  of  the  driving  function  G(t)  with  the  impulse  response  H(t).  That 
li» 

R(t)  -  H(t)*G(t)  ,  (H.4) 


where  *  denotes  convolution. 


What  is  usually  desired  is  to  obtain  tha  N  system  poles  and  residues 

of  (H.3).  This  must  be  done  by  deconvolution  of  the  driving  function 

(H.1)  from  (H.2).  If  the  r^  and  s^  are  known,  the  MfN  values  of  Sj,  contain 

the  M  values  s . .  The  g,  and  the  s.  are  known,  then  the  N  values  of  h  and 
J  J  j  k 


3k  can  be  determined  analytically. 


It  can  be  shown  that  the  response  function  can  be  written  in  terms  of 
the  driving  function  and  the  impulse  response  as 


M  N  \  ■  t  N  M  g.  s 

I  («i  I  rrr> • J  -  I  c-t  I  7 -£-> * 

J-l  3  k-i  s  k  k-1  k  >1  aJ  *k 


(H  .5) 


Hence,  the  response  function  residues  can  be  defined  as 

N  hk 

r,  -  8,  1  •  for  1  "  1  •  M 

1  J  k-1  8j  ak 


(H.6a) 


r,  -  -h  V  — — L-  ,  for  i  -  M+l  ,  N+M 

1  *  j-l  V  k 


(H.6b) 


From  (H.6b^  the  residues  of  the  impulse  response  are  written  as 


hk-|3: 

j-i  3j 


k  -  1  ,  N 


i  -  M+k 


(H.  7) 


Thus  the  res'  * tes  or  amplitudes  of  the  impulse  response  can  be  obtained  from 
the  known  va.  js  r^,  g^,  s^,  s^.  Of  course,  since  the  M  values  of  are 
known  and  the  M+N  values  of  s^  are  presumed  to  contain  the  M  values  of  , 
then  the  N  values  of  s^  can  be  determined  by  inspection. 


Prony's  method,  however,  will  not  necessarily  give  the  true  values 
of  the  driving  function  poles  in  the  extracted  poles  of  the  response 
function.  Hence  a  constrained  Prony’s  method  is  required  so  that  analyti¬ 
cal  deconvolution  of  (H.7)  cau  be  obtained. 

Work  performed  so  far  Implements  Method  1  which  is  described  below. 

It  was  found  that  this  method  works  perfectly  as  long  as  no  noise  is 
present  in  the  signal.  As  soon  as  noise  is  introduced  in  the  signal  and 
a  least-squares  method  is  used  the  matrix  (H.17a)  is  corrupted  with  noise 
which  through  the  least-squares  process  also  corrupts  the  constrained 
parameter.  Experiments  have  confirmed  this  observation.  This  suggests 
that  the  constraining  Method  1  may  work  well  in  noisy  data  if  the  curve¬ 
fitting  Prony's  method  is  used.  This  has  not  been  tested  as  yet. 

Method  2,  forcing  the  polynomial  root  solver  to  find  certain  roots, 
has  also  not  been  tested.  This  is  because  if  the  coefficients  are  all 
corrupted  by  noise  through  the  least-squares  procedure  then  subtracting 
any  given  poles  out  of  the  polynomial  will  force  the  remaining  poles  to 
carry  the  burden  of  all  the  noise. 

H.2  Method  1 


In  the  implementation  of  Prony's  method,  an  Nth  order  polynomial  is 
solved  for  its  roots.  The  order  of  the  polynomial,  N,  is  the  number  of 
poles  being  sought  in  the  transient  data.  If  the  coefficients  of  the 
polynomial  are  denoted  as  a  then  the  polynomial  can  be  expressed  as  per 
Reference  H.l  as 

aQ  ♦  ajZ1  +  a2Z2  +  ...  +  aNZN  -  0  (H.8) 


H-4 


where  is  usually  set  equal  Co  unity.  The  N  roots,  Z4  ,  are  defined  as 


with  s ^  being  the  poles  sought,  and  At  being  the  time  step  size  used  in  the 
analysis. 


If  the  value  of  one  or  more  of  the  poles  is  known  -  that  is,  we  know 
some  of  the  -  then  the  can  be  substituted  into  (H.B).  For  example, 
if  s^  or  Z^  is  known,  then  (H.B)  can  be  written  as 

Og  +  <*2^i  *  *  *  "**  aN^i  "  ®  •  (H10) 


The  N+l  polynomial  coefficients  a  are  solved  for  in  Prony's  method  by 
solution  of  the  difference  equation 


N-l 


p-0 


P  P+K 


-I 


N+K’ 


k-0,l,...,Y“l 


y-m-n 


(HJ1) 


The  and  I^(R  are  the  samples  of  the  transient  signal  being  analyzed  and 
M  is  the  total  number  of  samples  being  used.  The  value  of  M  must  be  at  least 
equal  to  2N  to  give  N  sets  of  equations  in  the  N  unknowns  ap.  However,  if 
the  value  of  a  pole  is  known,  then  one  of  the  N  aquations  can  be  equation 
H.10  and  N-l  equations  of  the  form  of  H.ll  can  be  written. 


If  L  poles  are  known  a  priori,  the  L  equations  of  the  form  of  (H.10) 
can  be  written  as 


N-l  „  „ 

V  a  zp«-z,  £  —  1,  2,  ,.,,L 

p  l  l 

p-0  v 


(H.12) 
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and  y»M-N-L  equations  can  be  written  in  the  form  of  (H.ll) 

N-l 

p-o  °p  Ip+K  "  "Ik+n  *  k"0,;L . Y_1  ’  Y"M”N“L  •  (h.13) 

Hence  there  are  still  y-M-N  total  equations  to  solve  for  the  N  values  of  a 
however  the  system  is  constrained  by  the  knowledge  of  the  location  of  L  poles. 
As  is  usually  done  in  Prony's  method,  if  M-2N  then  the  set  of  equations  is 
inverted  and  solved.  If  M>2N  then  a  pseudo-inverse  procedure  is  used. 


Using  the  matrix  notation  of  Reference  H.l,  if  M-2N  and  L  poles  are  known, 
then  we  solve  the  equation 


AB  -  C 


(H.  14) 


where  A  is  a  square  matrix  defined  as 

1  Z1 

1  Z- 


N-l-L 


N-L 


•  •  »  Z 

•  •  •  Z 


N-l 

1 

N-l 

2 


.N-l 


...  I 


N-l 


...  I 


N 


^-l+l  *•*  ^2N-L-2 


(H.  15a) 
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and  B  and  C  are  vectors  defined  as 


°N-1 


(H. 15b) 


(H. 15b) 


L2N-1-L 


If  M>2N  and  L  poles  are  known,  then  the  solution  takes  on  a  pseudo-inverse 
or  least-squares  form  as 


T  T 

A  AB  ■  A  C 


(H.  16) 


Where  A  la  the  transpose  of  the  matrix  A  and  A  la  now  a  rectangular  matrix 
of  the  form 

i  z:  z*  ...  zj’1 

l  z  z 2  2n-x 

i  i>2  u  2  •••  Z2 


M-N-l-L 


...  I\ 


...  I. 


...  I 


M-L-2 


(H. 17a) 


(H. 17b) 


Consider  the  simple  example  of  a  two  pole  system 

st  at 

f(t)  -  2e  1  +  3a  L 


where  s^  ■  0  and  s2  •  -A.  Assume  that  s^  Is  known,  thus  giving  Z ^ 
Also,  let  f(t)  be  sampled  at  At  -  0.2  seconds,  giving  a  data  set  as 


2.6057 


2.2722 


2.1223 


Using  the  constrained  Prony's  method  in  the  square  system  form  of  Equation 
(H.14)  and  (H.15)  yields 


or 


1 

1 

ao 

"l 

5 

3.3480 

-°1- 

_2.6057 

The  solution  of  this  sst  of  squations  givss  ■  1.4493  and  aQ  ■  0.4493  so 
that  tha  polynomial  can  ba  writtan  as 

Z2  -  1.4493  Z  +  0.4493  -  0  . 

Tha  roots  of  this  polynomial  ara  Z^  ■  1.0  and  Z 2  •  0.4493,  giving  polas  of 
8^^  •  0.0  and  s2  ■  4.0003.  Tha  arror  in  ia  dua  to  truncation  arror.  Nota 
that  tha  constralnad  pola  was  ratumad  axactly. 

H.3  Mathod  2 

Anothar  approach  which  can  ba  usad  to  constrain  Prony's  mathod  to  usa 
information  about  known  polas  is  tha  modification  of  tha  polynomial  root 
finding  routina  MULLER.  Onca  tha  unconstrainad  Prony's  mathod  calculatas 
tha  coafficiants  of  tha  polynomial  (H.8)  and  if  soma  of  tha  roots  of  tha 
polynomial  ara  known  a  priori,  than  tha  locations  of  thosa  roots  can  ba 
passad  to  MULLER  and  it  will  not  hava  to  saarch  for  thosa  roots.  Sines  tha 
root  finding  routina  la  vary  tlma  consuming,  tha  knowlsdgs  of  tha  location 
of  any  of  tha  roots  will  presumably  save  computation  time. 

A  possible  flaw  with  this  approach  is  that  MULLER  is  forced  to  pre¬ 
sume  roots  of  the  polynomial  when  those  exact  roots  may  not  ba  contained  in 
the  polynomial.  That  is,  if  the  polynomial  has  not  been  constrained  to 
contain  the  known  roots,  as  per  Section  H.2,  then  the  known  roots  will 


not  necessarily  be  contained  in  it.  Forcing  an  unconstrained  polynomial 
to  have  certain  roots  could  grossly  perturb  the  location  of  the  other  roots 
being  sought. 

H.A  Method  3 

The  obvious  solution  to  the  flaw  presented  in  Section  H. 3  is  to  use 
both  the  methods  of  H.2  and  H.3  simultaneously.  That  1st  the  polynomial 
is  constrained  to  contain  the  known  poles  as  outlined  in  Section  H.2.  Then 
the  polynomial  root  finding  routine  MULLER  can  be  modified  to  extract  the 
known  roots  from  the  polynomial  before  it  begins  its  search  for  the  unknown 
roots. 

It  is  felt  that  this  approach  will  give  the  best  accuracy  and  will 
spaed  up  the  calculations  since  the  order  of  the  polynomial  is  effectively 
reduced. 
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APPENDIX  I 

REVERSING  THE  WAVEFORM  IN  TIME  TO  ELIMINATE 
EXTRANEOUS  RESONANCES 


Reversing  the  waveform  in  time  dapanda  on  cha  atatiaticai  properties 
of  cha  noiaa  which  do  not  changa  whan  Cha  waveform  ia  ravaraad.  Tha  polaa 
chac  raaulc  from  cha  noiaa,  cherafora,  ara  not  alcarad  by  ravaraal  in  cima 
whlla  cha  crua  polaa  ara  nagacad  or  flippad  chrough  cha  origin. 

If  cha  noiaa  laval  ia  high  the  noiaa  polaa  only  approximately  remain 
cha  aame  under  Cima  ravaraal  and  tha  true  poles  only  approximately  reflect 
chrough  the  origin.  If  Cima  ravaraal  ia  attempted  for  curve-fitting  Prony's 
method  it  is  found  that  all  polaa  flip  precisely.  Therefore,  this  method 
ia  not  affective  for  curve-fitting  Prony's  method. 

(Curve-fitting  Prony'a  method  uaaa  tha  solution  to  cha  Inhomogeneous  system 
QXaq,  in  tha  notation  of  Volume  1,  where  Q  is  a  square  matrix.) 

Tha  waveform  of  Figure  1.1  waa  uaad  in  a  numerical  example  of  the  time-!* 
ravaraal  method.  Least-squares  Prony's  method  waa  applied  to  tha  waveform. 
Estimates  at  the  S-Plana  polaa  ware  found  using  the  inhomogeneous  solution 
which  is  defined  in  Volume  I,  Section  2  as 

\  "  [Q^r1  <fq. 

The  waveform  consists  of  100  samples  and  was  corrupted  with  uncorrelated, 
Gauaslan-diatributed  noise  with  a  standard  deviation  of  0.1.  Figure  1.2 
displays  the  poles  obtained  from  the  forward  and  reversed  waveforms.  The 
dimensions  of  the  (NXn) -dimensional  matrix  Q  in  this  example  ara  M*76 
and  r„"24.  Tha  estimates  of  the  true  poles  are  quite  accurate.  Figure  1.3 
displays  the  poles  obtained  for  tha  same  waveform  but  different  matrix 
dimensions.  The  matrix  dimensions  are  M*60  and  n*40.  In  this  case,  the 
faxtraneoua  poles  do  not  all  remain  in  the  left-half  plane  which  is 
attributed  to  using  an  overly  square  matrix. 
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If  the  matrix  Is  too  "thin",  of  M/n  »1,  the  estimates  of  the  true 
poles  become  Inaccurate.  If  the  matrix  is  too  "fat",  or  M/n  “  1,  the 
extraneous  poles  do  not  all  remain  in  the  left-half  plane.  The  matrix 
shape  where  M/n  «  .1  appears  to  be  the  best  shape  for  good  estimates  and 
keeping  the  extraneous  poles  in  the  left-half  plana. 
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Figure  X~3  Effect  of  an  overly-square  data  matrix. 


Therefore,  it  appears  that  if  we  choose  M/n  *  3,  then  time 
reversal  can  be  used  to  distinguish  between  the  true  poles  and  the 
extraneous  poles  in  relatively  high  noise  levels. 

The  above  two  examples  show  that  there  is  a  trend  for  the 
noise  poles  in  the  least-squares  Prony's  method  to  occupy  a  higher 
frequency  portion  of  the  S-plane  and  to  be  highly  damped.  This  trend 
was  studied  by  letting  the  least-squares  Prony's  method  operate  on  a 
waveform  consisting  of  only  Gaussian-dlstributed  uncorrelated  noise  - 
no  signal.  The  poles  resulting  from  this  example  are  highly  damped  and 
evenly  distributed  in  the  z-plane  approximately  around  a  circle  within 
the  unit  circle. 


This  behavior  can  be  explained  by  the  fact  that  the  polynomial 

coefficients,  except  which  is  set  to  one,  all  tend  to  zero  for  the 

least-squares  method  operating  on  uncorrelated  noise.  The  polynomial 

N 

then  tends  toward  z  ■  e  where  e  tends  to  zero.  The  roots  of  this 
polynomial,  z^,  have  magnitudes 


It  then  follows  that  the  poles  should  be  highly  damped,  since 
|z^|  +  0  if  | c |  *0,  and  that  the  poles  are  evenly  distributed  about  the 
z-plane.  However,  for  curve-fitting  Prony's  method  the  a^,  i  ■  0,  ..., 
N-l  do  not  ten  1  to  zero  and  indeed  this  phenomenon  is  not  observed  in 
tests  using  cuive- fitting  Prony’s  method  on  all  noise. 
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This  behavior  of  the  noise  poles  in  least-squares  Prony  seems  to 
remain  approximately  the  same  even  if  the  waveform  is  not  entirely  noise. 

The  trend  that  is  seen  is  that  the  noise  poles  occupy  the  higher  frequencies, 
are  highly  damped  and  evenly  distributed  between  the  higher  frequencies; 
while  the  true  poles  are  approximately  at  their  uncorrupted  locations 
and  occupy  the  lower  frequencies.  That  is,  it  appears  as  though  the  noise 
poles  are  "crowded"  away  from  the  lower  frequencies  or,  perhaps  more 
accurately,  the  lower  frequency  noise  poles  become  lower  frequency  signal 
poles  when  the  waveform  is  no  longer  entirely  noise. 


